1. Background

The HVT package is a collection of R functions to facilitate building topology preserving maps for rich multivariate data analysis, see Figure 1 as an example of a 3D torus map generated from the package. Tending towards a big data preponderance, a large number of rows. A collection of R functions for this typical workflow is organized below:

1. Data Compression: Vector quantization (VQ), HVQ (hierarchical vector quantization) using means or medians. This step compresses the rows (long data frame) using a compression objective.

2. Data Projection: Dimension projection of the compressed cells to 1D,2D or Interactive surface plot with the Sammons Non-linear Algorithm. This step creates topology preserving map (also called embeddings) coordinates into the desired output dimension.

3. Tessellation: Create cells required for object visualization using the Voronoi Tessellation method, package includes heatmap plots for hierarchical Voronoi tessellations (HVT). This step enables data insights, visualization, and interaction with the topology preserving map useful for semi-supervised tasks.

4. Scoring: Scoring new data sets or test data and recording their assignment using the map objects from the above steps, in a sequence of maps if required.

5. Temporal Analysis and Visualization: Collection of functions that analyzes time series data for its underlying patterns, calculation of transitioning probabilities and the visualizations for the flow of data over time.

6. Dynamic Forecasting: Simulate future states of dynamic systems using Monte Carlo simulations of Markov Chain (MSM), enabling ex-ante predictions for time-series data.

What’s New?

• The new update focuses on the integration of time series capabilities into the HVT package by extending its foundational operations to time series data which is emphasized in this vignette.

• The new functionalities are introduced to analyze underlying patterns and trends within the data, providing insights into its evolution over time and also offering the capability to analyze the movement of the data by calculating its transitioning probability and creates elegant plots and GIFs.

Below are the new functions and its brief descriptions:

• plotStateTransition: Provides the time series flowmap plot.
• getTransitionProbability: Provides a list of transition probabilities.
• reconcileTransitionProbability: Provides plots and tables for comparing transition probabilities calculated manually and from markovchain function.
• plotAnimatedFlowmap: Creates flowmaps for both self state and without self state scenarios.

2. Experimental setup

The Lorenz attractor is a three-dimensional figure that is generated by a set of differential equations that model a simple chaotic dynamic system of convective flow. Lorenz Attractor arises from a simplified set of equations that describe the behavior of a system involving three variables. These variables represent the state of the system at any given time and are typically denoted by (x, y, z). The equations are as follows:

\[ dx/dt = σ*(y-x) \] \[ dy/dt = x*(r -z)-y \] \[ dz/dt = x*y-β*z \] where dx/dt, dy/dt, and dz/dt represent the rates of change of x, y, and z respectively over time (t). σ, r and β are constant parameters of the system, with σ(σ = 10) controlling the rate of convection, r(r=28) controlling the difference in temperature between the convective and stable regions, and β(β = 8/3) representing the ratio of the width to the height of the convective layer. When these equations are plotted in three-dimensional space, they produce a chaotic trajectory that never repeats. The Lorenz attractor exhibits sensitive dependence on initial conditions, meaning even small differences in the initial conditions can lead to drastically different trajectories over time. This sensitivity to initial conditions is a defining characteristic of chaotic systems.

In this notebook, we will use the Lorenz Attractor Dataset. This dataset contains 200,000 (Two hundred thousand) observations and 5 columns. The dataset can be downloaded from here.

The dataset includes the following columns:

• X: The X-coordinate of a point in the Lorenz attractor. It represents the measure of strength of the circulation of fluid in the convective flow
• Y: The Y-coordinate of a point in the Lorenz attractor. It represents the measure of the temperature of the fluid at a given point in space, but it is the horizontal temperature distribution
• Z: The Z-coordinate of a point in the Lorenz attractor. It represents the measure of the temperature of the fluid at a given point in space, but it is the vertical temperature distribution
• U: U represents the velocity or speed of the system at a particular point in the attractor
• t: The time variable (discrete time steps which are incremented by approximately 0.0002500012500062 units of time) associated with each point in the Lorenz attractor. It indicates how much time has elapsed since the beginning of the simulation. Each value represents a specific point in time, and they are spaced apart by a fixed time interval.

3. Notebook Requirements

This chunk verifies the installation of all the necessary packages to successfully run this vignette, if not, installs them and attach all the packages in the session environment.

list.of.packages <- c("dplyr", "kableExtra", "plotly", "purrr", "data.table", "gridExtra", "grid", "reactable",
                    "reshape", "tidyr",  "stringr", "DT", "knitr", "feather","HVT","gganimate","gifski")

new.packages <-
  list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
if (length(new.packages))
  install.packages(new.packages, dependencies = TRUE, verbose = FALSE, repos='https://cloud.r-project.org/')
invisible(lapply(list.of.packages, library, character.only = TRUE))

4. Data Understanding

Here, we load the data. Let’s explore the Lorenz Attractor Dataset. For the sake of brevity, we are displaying only the first ten rows.

file_path <- ("./sample_dataset/lorenz_attractor.feather")
dataset <- read_feather(file_path) %>% as.data.frame()
dataset <- dataset %>% select(X,Y,Z,U,t)
dataset$t <- round(dataset$t, 5)
displayTable(head(dataset, 10))

Now, let’s try to visualize the Lorenz attractor (overlapping spirals) in 3D Space.

data_3d <- dataset[sample(1:nrow(dataset), 1000), ]
plot_ly(data_3d, x= ~X, y= ~Y, z = ~Z) %>% add_markers( marker = list(
                          size = 2,
                          symbol = "circle",
                          color = ~Z,
                          colorscale = "Bluered",
                          colorbar = (list(title = 'Z'))))

Now, let’s have a look at the structure of the Lorenz Attractor dataset.

str(dataset)

## 'data.frame':    200000 obs. of  5 variables:
##  $ X: num  0 0.0025 0.00499 0.00747 0.00995 ...
##  $ Y: num  1 1 1 0.999 0.999 ...
##  $ Z: num  20 20 20 20 19.9 ...
##  $ U: num  0 0.0005 0.001 0.0015 0.002 ...
##  $ t: num  0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 0.00225 ...

edaPlots(dataset)

edaPlots(dataset, output_type = 'histogram')

edaPlots(dataset, output_type = 'boxplot')

edaPlots(dataset, output_type = 'correlation')

edaPlots(dataset, time_column = "t", output_type = "timeseries")

5. Model Constructing and Visualization

The dataset is prepped and ready for constructing the HVT model, which is the first and most prominent step. Model Training involves applying Hierarchical Vector Quantization (HVQ) to iteratively compress and project data into a hierarchy of cells. The process uses a quantization error threshold to determine the number of cells and levels in the hierarchy. The compressed data is then projected onto a 2D space, and the resulting tessellation provides a visual representation of the data distribution, enabling insights into the underlying patterns and relationships.

Model Parameters

• Number of Cells at each Level = 100
• Maximum Depth = 1
• Quantization Error Threshold = 0.1
• Error Metric = Max
• Distance Metric = Manhattan / L1_Norm
• Dimension Reduction Metric = Sammon

NOTE: The compression takes place only for the X, Y and Z coordinates and not for U(velocity) and t(Timestamp). After training & Scoring, we merge back the U and t columns with the dataset.

hvt.results <- trainHVT(
  dataset[,-c(4:5)],
  n_cells = 100,
  depth = 1,
  quant.err = 0.1,
  normalize = TRUE,
  distance_metric = "L1_Norm",
  error_metric = "max",
  quant_method = "kmeans",
  dim_reduction_method = "sammon")

Let’s check out the compression summary.

summary(hvt.results)

NOTE: Based on the provided table, it’s evident that the ‘percentOfCellsBelowQuantizationErrorThreshold’ value is 0.02, indicating that only 2% compression has taken place for the specified number of cells, which is 100. Typically, we would continue increasing the number of cells until at least 80% compression occurs. However, in this vignette demonstration, we’re not doing so, because the plots generated from temporal analysis functions would become cluttered and complex, making explanations less clear.

Now, Let’s plot the Voronoi tessellation for 100 cells.

plotHVT(
  hvt.results,
  centroid.size = c(0.6),
  plot.type = '2Dhvt',
  cell_id = FALSE)

Figure 2: The Voronoi tessellation for layer 1 shown for the 100 cells in the dataset ’Lorenz attractor’

To understand how cell IDs are distributed across the map, we again plot Voronoi tessellation with cell_id = TRUE.

plotHVT(
  hvt.results,
  centroid.size = c(0.6),
  plot.type = '2Dhvt',
  cell_id = TRUE)

Figure 3: The Voronoi tessellation for layer 1 shown for the 100 cells in the dataset ’Lorenz attractor’ with Cell ID

6. Scoring

Once the model is constructed, the next step is scoring the data points using the trained HVT model. Scoring involves assigning each data point to a specific cell based on the trained model. This process helps map data points to their correct hierarchical cell without the need for forecasting.

scoring <- scoreHVT(dataset,
                    hvt.results,
                    child.level = 1)

Here we are displaying only the first 20 rows of the scored dataset.

summary(scoring)

These are the brief explanation of all the features in the above table.

• Segment Level: The tier or depth of a segment in the hierarchical structure.

• Segment Parent: The ID of the larger segment that this segment is part of in the level above.

• Segment Child: The IDs of smaller segments that are contained within this segment in the level below.

• n: The number of entities/data points from the uploaded dataset that are present inside that Segment child.

• Cell.ID: The ID of the child which is the result of sorted 1D Sammon’s.

• Quant.Error: The quantization error of that cell.

• centroidRadius: The maximum quantization error values as radius for anomalies.

• diff: The difference between centroidRadius and Quant.Error.

• anomalyFlag: The binary value that says the cell is an anomaly or not. (if the Quant.Error is greater than mad.threshold then it is = 1 (anomaly) else = 0 (not an anomaly))

Let’s look at the scored model summary from scoreHVT()

scoring$model_info$scored_model_summary

## $input_dataset
## [1] "200000 Rows & 5 Columns"
## 
## $scored_qe_range
## [1] "5e-04 to 0.3433"
## 
## $mad.threshold
## [1] 0.2
## 
## $no_of_anomaly_datapoints
## [1] 1706
## 
## $no_of_anomaly_cells
## [1] 18

The model info displays five attributes which are explained below:

1. Dimension of the dataset that is used in scoring.
2. The range of ‘Quantization Error’ after scoring of all the data points of the given dataset. (from minimum to maximum)
3. The value of Mean Absolute Deviation(MAD) Threshold.
4. The value of how many data points are anomalous. If the quantization error of a data point is greater than the ‘mad.threshold’, then it is considered as anomaly.
5. The value of how much cells have anomaly. If a cell has even 1 anomalous data point, the cell will be considered as anomaly_cell.

Now, let’s merge back the U and t columns from the dataset to the scoring output and prepare the dataset for ‘Temporal Analysis’ Functions.

temporal_data <- cbind(scoring$scoredPredictedData, dataset[,c(4,5)]) %>% select(Cell.ID,t)

7. Timeseries plot with State Transitions

The first novel function - plotStateTransition which is used to create a time series plotly object. The plot displays the movement of data points across the cells over time. Below is the function signature and its arguments.

plotStateTransition(
       df,
       sample_size,
       line_plot,
       cellid_column,
       time_column,
       v_intercept,
       time_periods)

• df - A data frame contains Cell ID and Timestamps.

• sample_size - A numeric value to specify the sampling value which ranges between 0.1 to 1. The highest value 1, outputs a plot with the entire dataset. Sampling of data takes place from the last to the first.

• line_plot - A Logical value. If TRUE, the output will be a timeseries plot with a line connecting the states according to the sample_size. If FALSE, a timeseries plot but without a line based on the sample_size will be the output.

• cellid_column - A character string specifying the column name of the Cell ID.

• time_column - A character string specifying the column name of the time stamp.

• v_intercept - A POSIXct/numeric object to draw a vertical line that denotes the ‘End of Training’ timestamp.

• time_periods - A list of vectors to draw the gray bars across the entire timeline. The bars represent the time periods with start and end timestamps.

plotStateTransition(df = temporal_data, 
                    cellid_column = "Cell.ID", 
                    time_column = "t",
                    sample_size = 0.2,
                    line_plot = TRUE)

For the demo of ‘sample_size’ argument, we are replicating the same plot above with the entire dataset.

plotStateTransition(df = temporal_data, 
                    cellid_column = "Cell.ID", 
                    time_column = "t",
                    sample_size = 1,
                    line_plot = TRUE)

8. Transition probability tables

The second novel function - getTransitionProbability which is used to calculate the transition probability of the current states to the next states including or excluding self-states. Below is the function signature and its arguments.

getTransitionProbability(
        df, 
        cellid_column, 
        time_column)
        type = "with_self_state"

• df - A data frame contains Cell ID and Timestamps.

• cellid_column - A character string specifying the column name of the Cell ID.

• time_column - A character string specifying the column name of the time stamp.

• type - A character string specifying how transition probability matrix should be calculated. Accepted values are ‘with_self_state’ or ‘without_self_state’. Default value is ‘with_self_state’.

trans_table <- getTransitionProbability(df = temporal_data,cellid_column = "Cell.ID",time_column = "t")
displayTable(trans_table, limit = 325)

• Current_State: The cell (out of 100 cells given in model training) in which the datapoint resides at a given time (t).

• Next_State: The cell (out of 100 cells given in model training) to which the datapoint moves at the next time unit (t+1).

• Relative_Frequency: The number of times that the datapoint moves from that Current_State to that Next_State.

• Transition_Probability: The probability calculated from the Relative_Frequency. Individual Relative_Frequency divided by the total of Relative_Frequency for a particular Current_State.

• Cumulative_Probability: The sum of the Transition_Probability for all the Next_State for a particular Current_State.

9. Reconciling transition probability using markovchain package

The third novel function - reconcileTransitionProbability which is used to reconcile the transition probability of the current states to the next states manually and markovchain function considering self states and without self states. Below is the function signature and its arguments.

reconcileTransitionProbability(
                df, 
                hmap_type = "All", 
                cellid_column, 
                time_column)

• df - A data frame contains Cell ID and Time stamps.

• hmap_type - A character string, if set to ‘without_self_state’, reconciliation plots for manual and markovchain for the highest transition probability excluding the self-state are given as output. If set to ‘self_state’, reconciliation plots for manual and markovchain for the highest transition probability considering the self-state is given as output and if set to ‘All’, plots including and excluding self-state are given as output.

• cellid_column - A character string specifying the column name of the Cell ID.

• time_column - A character string specifying the column name of the time stamp.

reconcile_plots <- reconcileTransitionProbability(df = temporal_data, 
                                                  hmap_type = "All", 
                                                  cellid_column = "Cell.ID",
                                                  time_column = "t")

Reconciliation plots of transition probability with self_state

The transition probability of one state staying in the same state is calculated using manual calculations and the markovchain function is plotted for comparison. The darker diagonal cells indicate higher probabilities of cells staying in the same state.

reconcile_plots[[1]]

Reconciliation table of transition probability with self-state

displayTable(reconcile_plots[[2]], limit = 217)

Reconciliation plots of transition probability without self-state

The transition probability of one state moving to the next state is calculated using manual calculations and the markovchain function and plotted for comparison. From all the next state transitions, the one with a higher probability is selected.

reconcile_plots[[3]]

Reconciliation table of transition probability without self-state

displayTable(reconcile_plots[[4]], limit = 217)

10. Flowmaps

The fourth novel function - plotAnimatedFlowmap which is used to create flowmaps and animations for the highest transition probability including and excluding self-states. Below is the function signature and its arguments.

plotAnimatedFlowmap(
         hvt_model_output, 
         transition_probability_df, 
         df, 
         flow_map = "All", 
         cellid_column, 
         time_column )

• hvt_model_output - The list object, which is the output from the trainHVT function.

• transition_probability_df - The probability list, which is the output from the getTransitionProbability function.

• df - A data frame contains Cell ID and Time stamps.

• flow_map - A character string. If set to ‘self_state’, the plot which shows the self-state by circles will be displayed. More the circle size, more probability of the cell stays in the same cell. If set to ‘without_self_state’, the plot which shows the next state by arrows will be displayed. The arrow head points the next cell to go from the cell which is pointed by the arrow tail. If set to ‘All’, two flow maps will be displayed. If set to NULL none will be displayed.

• cellid_column - A character string specifying the column name of the Cell ID.

• time_column - A character string specifying the column name of the Time stamp.

flowmap_plots <- plotAnimatedFlowmap(hvt_model_output = hvt.results,
                                     transition_probability_df =trans_table,
                                     df = temporal_data, 
                                     flow_map = 'All',
                                     cellid_column = "Cell.ID", 
                                     time_column = "t")

1. Flow map: Highest transition probability including self-state

The Circle size around the cell’s centroid represents self-state probability. More size, more probability of staying in the same cell.

flowmap_plots[[1]]

2. Flow map: Highest transition probability excluding self-states: Arrow size represents transition probability

The arrow size represents the Probability of the data to move to the next state. And the arrow directions point out to which cell it is moving next.

flowmap_plots[[2]]